On 05 April, 2024 I worked on(and if you are at the html knitted) this script.
The aim of the R Markdown script(or the output) is:
MPA Visualization - NS Germany
These files were downloaded from: https://sdi.eea.europa.eu/catalogue/srv/api/records/dae737fd-7ee1-4b0a-9eb7-1954eec00c65
Natura 2000 (vector) - version 2021 revision 1, Oct. 2022
🗺️ Zoom in to check the details in the OpenStreetMap in the background.
👆 Click on shapes to see German names.
Please follow the table of content (outline) for keeping
track of the steps. It must be on the left side.
DISCLAIMER:
Some code chunks could be eval=F to not to run the chunks but still
include them in reporting, or results=‘hide’might hide some results to
avoid cluttering.
Since the output document is shortened to demonstrate part of authors’
skill set.
Check the needed packages:
##print Rversion:#####
R.version.string
## [1] "R version 4.2.2 Patched (2022-11-10 r83330)"
#important for rmarkdown####
knitr::opts_chunk$set(echo = T, eval = T, fig.keep="all", cache = T) #DEFAULT:echo = T, eval = T, #eval=F does not run the chunk
knitr::opts_knit$set(root.dir = "~//")
#getwd()
#housekeeping
##packages####
options(rlib_downstream_check = FALSE) #dont check package downstream
#require(sp)
require(sf)
## Lade nötiges Paket: sf
## Linking to GEOS 3.11.1, GDAL 3.6.2, PROJ 9.1.1; sf_use_s2() is TRUE
require(rnaturalearth)
## Lade nötiges Paket: rnaturalearth
#library(knitr)
require(leaflet)
## Lade nötiges Paket: leaflet
#library(cowplot)
library(ggspatial)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.4
## ✔ forcats 1.0.0 ✔ stringr 1.5.0
## ✔ ggplot2 3.4.4 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.0
## ✔ purrr 1.0.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors
library(knitr)
library(kableExtra)
##
## Attache Paket: 'kableExtra'
##
## Das folgende Objekt ist maskiert 'package:dplyr':
##
## group_rows
library(units)
## udunits database from /usr/share/xml/udunits/udunits2.xml
##print sessionInfo:#####
#sessionInfo() #for reporting on package versions #commented out for privacy
Paths, parameters…
Loading files that were downloaded:
# Set the path to your GeoPackage file
gpkg_path <- paste0(shp_path,"EEA-Natura2000/eea_v_3035_100_k_natura2000_p_2021_v12_r01/GPKG/","Natura2000_end2021_rev1",".gpkg")
# Read the GeoPackage file
st_layers(gpkg_path)
## Driver: GPKG
## Available layers:
## layer_name geometry_type features fields
## 1 NaturaSite_polygon 27020 5
## 2 BIOREGION NA 28829 3
## 3 DESIGNATIONSTATUS NA 78751 5
## 4 DIRECTIVESPECIES NA 2315 9
## 5 HABITATS NA 152089 16
## 6 HABITATCLASS NA 139836 4
## 7 NATURA2000SITES NA 27027 23
## 8 OTHERSPECIES NA 334086 13
## 9 METADATA NA 4 2
## 10 IMPACT NA 200695 7
## 11 MANAGEMENT NA 32145 14
## 12 SPECIES NA 388303 19
## crs_name
## 1 ETRS89-extended / LAEA Europe
## 2 <NA>
## 3 <NA>
## 4 <NA>
## 5 <NA>
## 6 <NA>
## 7 <NA>
## 8 <NA>
## 9 <NA>
## 10 <NA>
## 11 <NA>
## 12 <NA>
n2kALL <- st_read(gpkg_path,layer = "NaturaSite_polygon")
## Reading layer `NaturaSite_polygon' from data source
## `/home/VTI_AD/oerey/03_Data/shapes/EEA-Natura2000/eea_v_3035_100_k_natura2000_p_2021_v12_r01/GPKG/Natura2000_end2021_rev1.gpkg'
## using driver `GPKG'
## Simple feature collection with 27020 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 746551.8 ymin: 905053.4 xmax: 6506127 ymax: 5298655
## Projected CRS: ETRS89-extended / LAEA Europe
#README<- st_read(gpkg_path,layer = "HABITATCLASS")
st_crs(n2kALL)
## Coordinate Reference System:
## User input: ETRS89-extended / LAEA Europe
## wkt:
## PROJCRS["ETRS89-extended / LAEA Europe",
## BASEGEOGCRS["ETRS89",
## ENSEMBLE["European Terrestrial Reference System 1989 ensemble",
## MEMBER["European Terrestrial Reference Frame 1989"],
## MEMBER["European Terrestrial Reference Frame 1990"],
## MEMBER["European Terrestrial Reference Frame 1991"],
## MEMBER["European Terrestrial Reference Frame 1992"],
## MEMBER["European Terrestrial Reference Frame 1993"],
## MEMBER["European Terrestrial Reference Frame 1994"],
## MEMBER["European Terrestrial Reference Frame 1996"],
## MEMBER["European Terrestrial Reference Frame 1997"],
## MEMBER["European Terrestrial Reference Frame 2000"],
## MEMBER["European Terrestrial Reference Frame 2005"],
## MEMBER["European Terrestrial Reference Frame 2014"],
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]],
## ENSEMBLEACCURACY[0.1]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4258]],
## CONVERSION["Europe Equal Area 2001",
## METHOD["Lambert Azimuthal Equal Area",
## ID["EPSG",9820]],
## PARAMETER["Latitude of natural origin",52,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",10,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["False easting",4321000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",3210000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["northing (Y)",north,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["easting (X)",east,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["Statistical analysis."],
## AREA["Europe - European Union (EU) countries and candidates. Europe - onshore and offshore: Albania; Andorra; Austria; Belgium; Bosnia and Herzegovina; Bulgaria; Croatia; Cyprus; Czechia; Denmark; Estonia; Faroe Islands; Finland; France; Germany; Gibraltar; Greece; Hungary; Iceland; Ireland; Italy; Kosovo; Latvia; Liechtenstein; Lithuania; Luxembourg; Malta; Monaco; Montenegro; Netherlands; North Macedonia; Norway including Svalbard and Jan Mayen; Poland; Portugal including Madeira and Azores; Romania; San Marino; Serbia; Slovakia; Slovenia; Spain including Canary Islands; Sweden; Switzerland; Türkiye (Turkey); United Kingdom (UK) including Channel Islands and Isle of Man; Vatican City State."],
## BBOX[24.6,-35.58,84.73,44.83]],
## ID["EPSG",3035]]
#transform projection for leflet #4326 is the EPSG code for WGS84
n2kALL<-st_transform(n2kALL, crs = 4326)
# Explore the structure of the spatial data
print(n2kALL)
## Simple feature collection with 27020 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -32.36806 ymin: 27.58768 xmax: 34.10149 ymax: 70.02586
## Geodetic CRS: WGS 84
## First 10 features:
## SITECODE SITENAME MS SITETYPE
## 1 PLH020044 Stawy Sobieszowskie PL B
## 2 PLH020088 Dalkowskie Jary PL B
## 3 PLH220055 Bunkier w Oliwie PL B
## 4 PLH020051 Irysowy Zagon koło Gromadzynia PL B
## 5 PLH020069 Las Pilczycki PL B
## 6 PLH020068 Muszkowicki Las Bukowy PL B
## 7 PLH020070 Sztolnia w Młotach PL B
## 8 PLH020078 Kumaki Dobrej PL B
## 9 PLH020084 Dolina Dolnej Baryczy PL B
## 10 PLH020087 Gałuszki w Chocianowie PL B
## INSPIRE_ID geom
## 1 PL.ZIPOP.1393.N2K.PLH020044 MULTIPOLYGON (((15.67494 50...
## 2 PL.ZIPOP.1393.N2K.PLH020088 MULTIPOLYGON (((15.87839 51...
## 3 PL.ZIPOP.1393.N2K.PLH220055 MULTIPOLYGON (((18.55123 54...
## 4 PL.ZIPOP.1393.N2K.PLH020051 MULTIPOLYGON (((16.35672 51...
## 5 PL.ZIPOP.1393.N2K.PLH020069 MULTIPOLYGON (((16.94962 51...
## 6 PL.ZIPOP.1393.N2K.PLH020068 MULTIPOLYGON (((16.95908 50...
## 7 PL.ZIPOP.1393.N2K.PLH020070 MULTIPOLYGON (((16.54941 50...
## 8 PL.ZIPOP.1393.N2K.PLH020078 MULTIPOLYGON (((17.21745 51...
## 9 PL.ZIPOP.1393.N2K.PLH020084 MULTIPOLYGON (((16.5485 51....
## 10 PL.ZIPOP.1393.N2K.PLH020087 MULTIPOLYGON (((15.71207 51...
# View basic summary information
summary(n2kALL)
## SITECODE SITENAME MS SITETYPE
## Length:27020 Length:27020 Length:27020 Length:27020
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
## INSPIRE_ID geom
## Length:27020 MULTIPOLYGON :27020
## Class :character epsg:4326 : 0
## Mode :character +proj=long...: 0
List of member states in n2kALL:
unique(n2kALL$MS)
## [1] "PL" "EE" "MT" "SK" "DE" "PT" "DK" "SI" "RO" "SE" "IE" "HR" "FI" "CZ" "ES"
## [16] "FR" "IT" "LV" "AT" "LT" "CY" "LU" "NL" "BE" "BG" "GR" "HU"
#unique(n2kALL$SITETYPE)
#filterout only DE
n2kde <- n2kALL[n2kALL$MS == "DE", ]
#class(n2kde)
#plot what we have (takes time!)
#plot(n2kde)
rm(n2kALL)
cat("Total number of polygons after filtering:",length(st_geometry(n2kde)))
## Total number of polygons after filtering: 5200
# Define the bounding box based on your specified range
bbox <- st_bbox(c(xmin = mgmt_lon[1], xmax = mgmt_lon[2], ymin = mgmt_lat[1], ymax = mgmt_lat[2]), crs = 4326)
# Convert the bounding box to a simple features (sf) polygon
bbox_polygon <- st_as_sfc(bbox)
#plot(bbox_polygon)
leaflet() %>%
addTiles() %>%
addPolygons(data= bbox_polygon)
# Set the path to your GeoPackage file
DEeez_path <- paste0(shp_path,"EEA-Natura2000/DEeez/")
DEeez_sf <- st_read(DEeez_path)
# Crop 'coast_sf' to the specified bounding box
DEeez_sf <- st_intersection(DEeez_sf, bbox_polygon)
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
leaflet() %>%
addTiles() %>%
addPolygons(data= DEeez_sf, color = "black") #%>%
#addPolygons(data= DEeez_sf_buffered, color = "purple")
# Crop 'coast_sf' to the specified bounding box
n2kdeNS <- st_intersection(n2kde, bbox_polygon)
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
leaflet() %>%
addTiles() %>%
addPolygons(data= n2kdeNS, color = "black")
rm(n2kde)
print(paste("Total number of polygons after filtering:",length(st_geometry(n2kdeNS))))
## [1] "Total number of polygons after filtering: 208"
# Crop 'coast_sf' to the DEeez (but it also merges)
n2kdeEEZ <- st_intersection(n2kdeNS, DEeez_sf)
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
leaflet() %>%
addTiles() %>%
addPolygons(data= n2kdeEEZ, color = "pink")
rm(n2kdeNS)
print(paste("Total number of polygons after filtering:",length(st_geometry(n2kdeEEZ))))
## [1] "Total number of polygons after filtering: 43"
We want to focus on only the biggest of MPAs in german EEZ and territorial waters. I decided to calculate the area size and sort them accordingly.
#add column with area
n2kdeEEZ$area <- st_area(n2kdeEEZ) #Take care of units
class(n2kdeEEZ$area)
## [1] "units"
#check unit
attr(n2kdeEEZ$area, "units")
## $numerator
## [1] "m" "m"
##
## $denominator
## character(0)
##
## attr(,"class")
## [1] "symbolic_units"
#it is meter squares
#turn it into kmsquare by /1000000
n2kdeEEZ$area2 <- as.numeric(n2kdeEEZ$area/1000000)
summary(n2kdeEEZ$area2)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 0.160 9.902 640.751 109.824 5289.394
n2kdeEEZ %>%
ggplot(aes(x = area2)) +
geom_histogram(binwidth = 10, fill = "black", color = NA) +
labs(x = "Area", y = "Frequency", title = "Histogram of Area")+
theme_minimal()
knitr::kable(as.data.frame(n2kdeEEZ) %>%
summarize(
TotalRows = n(),
RowsWithAreaGreaterThan100 = sum(area2 > 100),
RowsWithAreaGreaterThan50 = sum(area2 > 50),
RowsWithAreaGreaterThan40 = sum(area2 > 40),
RowsWithAreaGreaterThan20 = sum(area2 > 20),
RowsWithAreaGreaterThan10 = sum(area2 > 10)
))
| TotalRows | RowsWithAreaGreaterThan100 | RowsWithAreaGreaterThan50 | RowsWithAreaGreaterThan40 | RowsWithAreaGreaterThan20 | RowsWithAreaGreaterThan10 |
|---|---|---|---|---|---|
| 43 | 11 | 15 | 15 | 17 | 21 |
leaflet() %>%
addTiles() %>%
addPolygons(data= n2kdeEEZ, color="black")
smallestareas<- n2kdeEEZ %>% dplyr::filter(area2<50)
leaflet() %>%
addTiles() %>%
addPolygons(data= smallestareas, color="black",
popup = ~SITENAME)
These must be partial shapes due to cutting out of the DE-EEZ and territorial waters.
n2kdeEEZ %>% nrow()
n2kdeEEZ %>% dplyr::filter(area2>50) %>% nrow()
biggestareas<- n2kdeEEZ %>% dplyr::filter(area2>50)
leaflet() %>%
addTiles() %>%
addPolygons(data= biggestareas, color="black")
I decided to show only the shapes that are bigger than 50 kilometer square.
cat("Final number of polygons after filtering:",length(st_geometry(biggestareas)))
## Final number of polygons after filtering: 15
Since types overlap. Click to read the name. Use the legend on top right to choose site type.
unique(n2kdeEEZ$SITETYPE)
## [1] "B" "C" "A"
# split data
biggestareas.split <- split(biggestareas, biggestareas$SITETYPE)
# Define the color palette
pal <- colorFactor("Set1", unique(biggestareas$SITETYPE))
# Initialize Leaflet map
l <- leaflet() %>% addTiles()
# Iterate over each group in biggestareas.split
names(biggestareas.split) %>%
purrr::walk(function(group) {
l <<- l %>%
addPolygons(data = biggestareas.split[[group]],
color = ~pal(biggestareas.split[[group]]$SITETYPE),
fillOpacity = 0.5,
fillColor = ~pal(biggestareas.split[[group]]$SITETYPE),
popup = ~SITENAME,
group = group)
})
# Add layers control
l %>%
addLayersControl(
overlayGroups = names(biggestareas.split),
options = layersControlOptions(collapsed = FALSE)
)
n2KdeEEZbiggerthan50km2<-paste0(data_out,"/n2KdeEEZbiggerthan50km2")
# Write the sf object to a shapefile
st_write(biggestareas, n2KdeEEZbiggerthan50km2, driver = "ESRI Shapefile")
knitr::kable(as.data.frame(biggestareas) %>% dplyr::select(c("SITENAME","SITETYPE","SITECODE","area2")),
caption = "FINAL DE MPA SHAPES") %>%
kable_styling(full_width = FALSE)
| SITENAME | SITETYPE | SITECODE | area2 |
|---|---|---|---|
| Steingrund | B | DE1714391 | 173.72694 |
| Borkum-Riffgrund | B | DE2104301 | 618.11929 |
| Doggerbank | B | DE1003301 | 1676.77056 |
| Nationalpark Niedersächsisches Wattenmeer | B | DE2306301 | 2574.97468 |
| Sylter Außenriff | B | DE1209301 | 5289.39409 |
| NTP S-H Wattenmeer und angrenzende Küstengebiete | B | DE0916391 | 4283.49474 |
| Helgoland mit Helgoländer Felssockel | B | DE1813391 | 54.11029 |
| Unterelbe | B | DE2018331 | 78.60430 |
| Hamburgisches Wattenmeer | C | DE2016301 | 133.62144 |
| Unterems und Außenems | B | DE2507331 | 51.94993 |
| Schleswig-Holsteinisches Elbästuar und angrenzende Flächen | B | DE2323392 | 86.02672 |
| Ramsar-Gebiet S-H Wattenmeer und angrenzende Küstengebiete | A | DE0916491 | 4305.92962 |
| SPA Östliche Deutsche Bucht | A | DE1011401 | 3120.12815 |
| Seevogelschutzgebiet Helgoland | A | DE1813491 | 1605.77646 |
| Niedersächsisches Wattenmeer und angrenzendes Küstenmeer | A | DE2210401 | 3351.12014 |
This is just a plotting exercise. Please contact me for details: serra.oerey@uol.de 🧐
End of the document. - Serra Örey 💜